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Abstract 

We describe in detail how perturbations due to the planets can cause a sub- 
population of WIMPs captured by scattering in surface layers of the Sun to 
evolve to have orbits that no longer intersect the Sun. We argue that such 
WIMPs, if their orbit has a semi-major axis less than 1/2 of Jupiter's, can 
persist in the solar system for cosmological timescales. This leads to a new, 
previously unanticipated WIMP population intersecting the Earth's orbit. 
The WIMP-nucleon cross sections required for this population to be significant 
are precisely those in the range predicted for SUSY dark matter, lying near 
the present limits obtained by direct underground dark matter searches using 
cyrogenic detectors. Thus, if a WIMP signal is observed in the next generation 
of detectors, a potentially measurable signal due to this new population must 
exist. This signal, lying in the keV range for Germanium detectors, would 
be complementary to that of galactic halo WIMPs. A comparison of event 
rates, anisotropics, and annual modulations would not only yield additional 
confirmation that any claimed signal is indeed WIMP-based, but would also 
allow one to gain information on the nature of the underlying dark matter 
model. 



I. INTRODUCTION 

It is by now firmly established that the dynamics of galaxies and clusters of galaxies is 
governed by the existence of large amounts of dark matter, with an average density sufficient 
to result in a total mass density in excess of 10% of the closure density today |1]]. At the 
same time, constraints from Big Bang Nucleosynthesis suggest that the abundance 

of baryons, including dark baryons, is not likely to be sufficient to account for all of this 
material. Thus, one is led to the possibility of non-baryonic dark matter, composed of a 
dissipationless gas of very weakly interacting elementary particles. The data from structure 
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formation also suggests that most of the dark matter must be "cold" , namely non-relativistic 
at the time fluctuations on the scale of galaxies first could begin to condense. 

For these reasons, a favored candidate for dark matter is a so-called WIMP, Weakly 
Interacting Massive Particle. The best motivated among the possibilities involves super- 
symmetric extensions of the standard model. In these models, the lightest supersymmet- 
ric partner of ordinary particles, usually the neutralino, can be absolutely stable. More- 
over, in order to resolve various naturalness problems in the Standard Model, the mass 
scale of the neutralino is expected to be comparable to the weak symmetry breaking 
scale. As a result, there is a dynamical argument that naturally leads to a remnant neu- 
tralino WIMP abundance comparable to the closure density today, (e.g. see 0): The 
fraction of the closure density provided by cold relics left over after out-of-equilibrium 
annihiliation of an initially thermal distribution of particles, say X, is on the order of 

^x/Pciosure = ~ 10"^''' cm^/i~^/((Ta(f /c)) whcre (Ta is the annihilation cross section 
(and where Qx = Px/Pciosure, and h = i^o/lOO km/s/Mpc) |^J^.As a result, the typical 
annihilation cross section needed for leaving a density of massive relics comparable to the 
closure density is a weak scale cross section: q;^(100 GeV)~^ ~ (a/lO"^)^ x 4 x 10"'^^ cm^. In 
this paper, we shall assume that the dark matter is indeed made of such weakly interacting 
massive particles (WIMPs) . Though most of our work depends only on the general assump- 
tion of WIMPS (with mass 10 GeV < mx < 1000 GeV), we shall give numerical estimates 
for the effects we discuss by sampling the parameter space of the Minimal Supersymmetric 
Standard Model (MSSM) in its various forms 0. 

Several searches for these WIMPs are underway, using either direct (energy deposition in 
laboratory samples) or indirect (e.g. looking for products of WIMP annihilation) techniques. 
Direct searches look for recoil events (with associated phenomena such as heat deposition, 
ionisation, etc.) due to the elastic scattering of WIMPs on nuclei. These searches present 
many experimental challenges because: (i) the rate of expected signals is small (less than a 
few events/day/kg), (ii) the recoil energies are also small (typically from 1 to 30 keV), and 
(iii) the background rate is comparatively very high. 

Any a priori theoretical information on the expected signals is therefore very important 
for separating true signals from background events. In the present paper, we shall describe 
in detail the derivation and characteristics of our previously announced result that a 
heretofore unexplored aspect of the dynamical history of WIMPs in the solar system may lead 
to a new population of WIMPs whose signals lead to a peculiar signature in the differential 
energy spectrum of recoil events, and, probably, peculiar anisotropy and annual modulation 
features. The new signals discussed here correspond to an excess of recoil events in an energy 
window on the keV energy scale. The ratio of the rate of these events (per day, per kg and per 
keV) to the standardly expected one is proportional to an average scattering cross section of 
the WIMP, weighted over elements in the Sun. Sampling over SUSY parameter space, we find 
that this event rate ratio can be larger than a factor of 2. Most important perhaps, we find 
that it is for cross sections that lie just below current experimental limits, namely for those 
WIMPs that are the prime target of the next generation of underground detectors, that the 
new signal we describe is maximal. Hence, if WIMPs are discovered in the next generation of 
detectors, the signal we propose may be one of the best ways of demonstrating its origin, and 
probing the characteristics of the WIMP dark matter responsible for it. Finally, another 
signal associated to the new dynamical class of WIMPs that we discuss here might be a 
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significant increase in the (indirect) neutrino signal expected from WIMP anniliilations in 
tlie Eartli. Indeed, limits from such searches might possibly further constrain the region of 
SUSY dark matter parameters that remain viable. 

The outline of this paper is as follows: First, we derive in section II the differential 
capture rate by the Sun for WIMPs that might eventually form part of the population of 
interest here. Next (section III), we derive in detail the dynamical equations that govern the 
possible diffusion of this population into bound solar system orbits that no longer intersect 
the Sun. In the next section, we combine the results of the two preceding discussions to 
estimate the capture rate of long- surviving solar-system bound WIMPs. Based on this, we 
can then estimate (in section V) the present phase-space distribution of such WIMPS in 
the region of the Earth. We are then able, in section VI, to explore the possible observable 
signals from this new population, in terms of the differential event rate per kg per day per 
keV of scattering on target nuclei. In section VII we explore these results in the context 
of realistic SUSY WIMPS. By sampling the allowed SUSY phase space, using accelerator 
constraints, and using the "Neutdriver" code to determine remnant WIMP densities and 
calculate capture and scattering cross sections on various targets, we demonstrate that the 
direct scattering signal from this new population should in principle be detectable if WIMPs 
lie near the current bounds. Finally, in the concluding section we examine other possible 
signatures for this distribution, and discuss future work that would be useful to perform in 
light of the results we present here. 



II. DIFFERENTIAL CAPTURE OF WIMPS BY THE SUN 

Direct searches of WIMPs assume that the recoil events are due to the scattering on 
a nucleus of a WIMP coming directly from a "standard" galactic halo, with a local mass 
density (around the solar system) px ~ 0.3GeV/cm^, and a roughly Maxwellian velocity 
distribution (in the galactic inertial frame) characterized by v^ras ~ 270 km/sec. In this 
paper, we shall study a particular class of WIMPs that underwent the following dynamical 
history: 

(i) coming from the galactic halo, they are scattered by nuclei in the outskirts of the Sun 
into very elliptic, bound orbits with semi-major axis a ~ 1 astronomical unit (AU); 

(ii) being perturbed by gravitational interaction with the planets, they diffuse out of the 
Sun and stay bound in the inner solar system during its entire age, ~ 4.5 Gyr; 

(iii) finally, they form a class of low velocity WIMPs, with typical barycentric velocities 
vx ~ "yEarth ~ 30km/s, which can either undergo a scattering event with a nucleus in a 
dark-matter detector, or be scattered by a nucleus in the Earth to be ultimately accreted to 
finally annihilate at the Earth's center. 

In this Section, we consider the first step of this process: the scattering by nuclei in 
the Sun into an AU-scale bound orbit. The essential new feature with respect to previous 



analyses of capture of WIMPs by the Sun pHlO[ will be to concentrate on WIMPs that 
graze the Sun, and lose just enough energy to stay in Earth-crossing orbits. To estimate 
the number of WIMPs susceptible to be sufficiently perturbed by the small gravitational 
interaction of planets, we need first to derive the differential capture rate, per energy and 
per angular momentum, of WIMPs by the Sun. We shall ultimately be interested only in the 
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small fraction of WIMPs which have angular momenta in a small range Jmin ^ J ^ Js where 
Js is the angular momentum for a WIMP exactly grazing the Sun. The result we need is a 
simple generalization of results previously derived in the literature 0,0. However, for the 
benefit of the general reader we shall present a self-contained derivation starting essentially 
from scratch. 
Let 

dnx = /oo(Voo) Xoo V^o = /loc(xioc, Vioc) d^ Xioc d-^ Vioc , (2.1) 

be the phase-space distribution of WIMPs, et infinity (in the galactic halo), and locally 
within the solar system. We shall always refer WIMP velocities to a frame attached to the 
Sun. Let denote the (vectorial) velocity of the Sun with respect to a galactic frame. 
The galactic WIMP velocity would be v^^^^jj^p = vwimp + ^s- The distribution function at 
infinity is taken to be Maxwellian, 

/-(^-) = ^^37^ \ -2 j (2-2) 

where = ffrms = | (^^oo) is an rms "planar" velocity. As the class of WIMPs studied here 
does not depend much on any eventual cut-off of foo beyond some "evaporation" velocity, we 
shall work with the simple exponential form (|2.2| ). In the body of the text we shall assume 
the following standard values for the parameters px = nx mx, Vo, and vs- 

^|andard _ _ g.S GeV cm'^ , y^t^-dard _ 220 kms"' , Vs = 220 kms-^ . (2.3) 

At the end, we shall comment on the effect of changes around the standard values for px 
and Vo- 

Liouville's theorem tells us that /ioc(xioc, vioc) is constant along the free motion of 
WIMPs. Therefore, at any point along an incoming trajectory 

/loc (Xloc, Vioc) foo [Vqo (Xloc, Vioc)] , (2.4) 

where Voo(xioc, vioc) is the incoming velocity at infinity of the WIMP observed locally, at 
position Xloc, with velocity vioc, within the solar system (before it undergoes any scatter- 
ing event). In spherically symmetric problems, only the angular average /loc (xioc, vioc) = 
/oo [foo (xioc, Vioc)] matters, with 
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47r^/^ VoVs Voc 



(2.5) 



It is standard to write the differential scattering cross-section of the WIMP X onto the 
nucleus of atomic number A as 

cia^ = aAFj(Q)^, (2.6) 

where Q = -^before — -^'after IS the energy transferred during the scattering, Fa{Q) is a form 
factor, and dQcm = sin ^cm d9cm dfcm is the scattering solid angle element in the center of 
mass frame. We shall take an exponential form factor 
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F\Q) = ew{-Q/QA), 



(2.7) 



where the nucleus-dependent quantity Qa will be discussed below. Note that cr^ (often 
denoted cx^) in Eq. ( p.6| ) is by definition independent of the scattering angle. 

Let us consider a volume element (i^x in the Sun, containing r;,^(x)(i^x nuclei of 
atomic number A. The differential flux of WIMPs impinging on this lump of scatterers 
is (i^vioc /ioc(x, vioc) |vioc|. Therefore, the corresponding differential number per second of 
scattering events (within the center-of-mass scattering solid angle dVLcm) reads 

dNA = (i^xnA(x) d\ioc /ioc(x, vioc) |vioc| cta F\{Q) . (2.8) 

47r 

We wish to sort out this scattering rate according to the distribution in outgoing semi-major 
axis a and (specific) angular momentum J. The conserved energy (kinetic plus potential) 
of the WIMP before the scattering event is 

^before = ^ ^xi^l^^ - t^esc(^)) = \^X^lo, (2-9) 

where vioc = Vbefore denotes the local velocity before the collision, and where fesc(^) = 
2U{r) = +2 J d^x.' GArp(x')/|x — x'l is the escape velocity at the radius r within the Sun {G^ 
denoting the Newtonian gravitational constant). The conserved energy after the collision 
reads 

F 1 ^2 2 f ^^ GxmxMQ _ 1 

Rafter = ^ "^X (v^f^ " ^csc(^)) = ^ = --TUx a , (2.10) 

where Mq denotes the mass of the Sun, and where a is a shorthand for G^MQ/a. By 
the standard laws of nonrelativistic elastic collisions (see, e.g., '^^), the velocity after the 
collision, and therefore the energy transfer, are linked to the cm. scattering angle 6cm by 



after loc 



l-^/?^(l-COS^c 



TP T7< ^ pA 2 ^ COs6'cm 
, — -^before — -Chaffer — - '^X P+ ^loc ^ 5 



(2.11) 

where we define (following Ref. [Q) 



i^i^ , ? (2.12) 
[mx ± vtlaY 

Note that the maximum value of is one, which is reached when the mass of the WIMP 
matches the mass of the nucleus: mx = mA- For a given incoming local velocity vioc, 
Eqs. dH), O give 



Q = l ^x{vl^ - vlJr) + a) = ^ mxivl + a) . (2.13) 

On the other hand, Eq. (|2.11|) relates the energy transfer Q to the cm. scattering solid 
angle. Finally, we have 
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d^cm _ , / 1 - COS gem A _ dQ _ da 



loc 



where we recall that a = Gn Mq/q. 

Inserting Eq. ( p.l4| ) in Eq. (|2^ ) yields 



diV^ = rf^xn^(x) ^'^'°- ^--(^'^'°-) ir2(g) rf« , (2.15) 



'loc 



where the last factor is a step function (0(x) = 1 if a; > 0, and vanishes if x < 0) taking 
care of the inequality on vioc or t>^(x, vioc) = ffoc — ^esc(^) entailed by the constraint (1 — 
cosecm)/2 < 1 , i.e. 2Q/mx = v^^ + a < iS^vl^ = (3^{vl^ + vl^{r)). Using the identity 
(1//?+) — = 1, easily checked from the definitions ( p.l2| ), this constraint leads to the 

step function 




(2.16) 



The phase-space distribution of the WIMPs at infinity, Eq. ( p.2|) , is anisotropic. If 
the overall capture mechanism discussed in this paper were spherically symmetric (i.e. a 
capture by a spherical Sun followed by a spherically symmetric exit mechanism), it would be 
exact to replace, for purposes of capture calculations, the (anisotropic) incoming phase-space 
distribution, Eq. ( p.2|) , by its (isotropic) angular average, Eq. ( |2.5| ). Actually, the overall 
capture mechanism discussed here is not spherically symmetric (even in the (excellent) 
approximation of an exactly spherical Sun), because, as we shall see in the next section, 
two angular parameters {i and g), related to the spatial orientations of Vafter and x with 
respect to the ecliptic plane, modulate the efficiency of extraction of the WIMPs captured 
by the Sun. However, we expect that, because of the partial randomisation of the incoming 
directional quantities by scattering on a spherical Sun, the effects linked to the overall lack 
of spherical symmetry, i.e. the correlation effects between the incoming anisotropy (linked to 
the direction of V5) and the outgoing one (linked to the ecliptic plane) are small. Therefore, 
we shall henceforth neglect them, for purposes of capture calculations, and replace the 
anisotropic distribution, Eq. (|2^ ), by its angular average, and, correspondingly, the local 
distribution /ioc(x, vioc) by its angular average f^oc which, thanks to Liouville's theorem 
( p.4|) , is simply equal to [voo{r,v\oc)]- We can then make use of the spherical symmetry 
of /iqj, and n^(x) = nA{r) to simplify the problem of sorting out the differential rate ( p. 15 ) 



according to the distribution in the outgoing (specific) angular momentum J = |x x Vafterl- 
(We thank A. Gould for suggesting this simplification). Let 6 denote the colatitude of 
Vafter with respect to the radial direction taken as z-axis, i.e. J = rVafter sin 9. Since an 
isotropic distribution is uniform in cos 6', and J oc sin 6*, the distribution in is given by 
the normalized measure 

dJ'{2JlJ-\l - jyjlJ-'^'Qj, (2.17) 

where 

^max(r, a) = r(fgg^(r, a) - af^^ = rt^after , 6j = 6(^max(r, a) - J) . (2.18) 
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Here, the theta function Bj takes care of the constraint sin 6' < 1. This leads to a differential 
scattering rate of the form 



A-KViocdVloc /loc(^7 ^^Ic 



'1 



J2 



^^/^F](g)e„ejc/arfj^ (2.19) 



Performing the integral over J in Eq. (|2.19|) over the range 



dNA\ 



j>.j„ 



d y^UAir) OA 1 



(2.20) 



in which Q should be replaced by its expression in terms of fj^^, or w^, Eq. ( |2.13| ). Here, 
® Jmin = ['^ "f^after " <^min] = © [<^max('", " <^mm] gives a lower bouud to the admissible values 
of r = |x|. The integrals over (i^x = Anr'^dr and dvioc are uncoupled. Using the exponential 
form factor (E.TI) we then define the function 



KA{r, a) = — —J / Atcvioc dvioc /loc exp 

= — -T4- / 4:71 dvoo /oo(^^oo) exp 
In the second form, simplified by taking Voo{r, vioc) = ( 



mx 



2Qy 



i^L - vLir) + a) 



Or 



mx_ 
'2Qa 



(2.21) 



-'loc 



we used f loc (if loc 



3c(^))^ as integration variable, 
dvoo and Liouville's theorem. Note that the step function Bq, limits 
to the range < < l3^{vi^{r) - a/ (3^). With the definition (^i^, the result 
can be finally written as (after integration over the Sun) 



dNA 



da 



nx 

Vo 



d^i^nAir) OA 1 



/2 



1/2 



r>r„ 



i^^(r, a) 



(2.22) 



Evidently, the sum over all the (significant) values of A (atomic number) present in the 
Sun must be ultimately performed. Here, the minimum radius rmin (perhilion distance) is 
defined in terms of the minimum angular momentum Jmin by r^^^ {vl^Sj-aa^ ~ «)^''^ = </min- 
This formula yields the result needed for our purpose. Namely, the rate with which WIMPs 
scatter on nuclei with atomic number A to end up into bound solar orbits with semi-major 
axis within [a, a + da\ (corresponding to [a, a + da\ with a = Gfq Mq/o), and with specific 
angular momentum J > Jmin- 

The limit Jmin — ^ 0^ would collect all the WIMPs scattered at any point within the Sun, 
even with very small perihelion distances rmin, i-e. passing very near the center of the Sun. 
In our subsequent work we shall be interested in the opposite limit Jmin ^ {Rs f^esdRs))' 
corresponding to orbits which graze the Sun (radius Rs) and barely penetrate it. The usual 
total capture rate is obtained from Eq. (|2.22|) by setting Jmin = and by integrating over 
a = GMq/q. We shall be interested in values Jmin — Rs'fesdRs) and a ~ lAU, i.e. 

v"^ where ve = 29.8 km/s is the Earth orbital velocity. Note that 
■^esc ^ « SO that we can, with a very good approximation, neglect 



a ~ GM0/(1AU) - 
for such values of a, 
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a with respect to v\^^ both in the square-root factor in Eq. (|2.22|) and in the definition of 
'^min: '^min "yescl^min) — <^min- Further, for Suu-graziug orbits we can approximate the radial 
dependence of the escape velocity by: fesc(^) — '^GMq/t. Within these approximations, the 
differential rate [(iA^^/rfa] j^.^ reads 



dNA 



da 



— / d^xnA(r) a A { 1 

Vo J r>r min 



1/2 



(2.23) 



Note also that for the values of a we are interested in, typically a ~ f|; <^ f^, so that 
the function KA{r,a) is nearly independent of a. Let us finally write more explicitly the 
result ( p.21|) . Inserting the explicit (angle-averaged) shifted Maxwellian spectrum (p.5| ) into 
Eq. ( p.21|) , one can express Ka in terms of the error function 



X{x) = dye 



TC 



1/2 



erf (x) 



(2.24) 



To do this one must use as integration variable x = v^l + Vao/vo (where a a is defined 
below). Setting 



mx vl 



Vs 



2Q, 



(2.25) 



and, with A{r) > 0, 



a 



\ p. 



A ' 



(2.26) 



yields 



KA{r, a) 



exp (-f^) expi-r^jaA] 



X (Air) -r]a)-X (^(r) +Va) + 2x (Va) ■ (2.27) 



TTy^/3^il+dA)Va 

In the applications below we shall use the standard value for the coherence energy Qa 



= 10"^^ cm 



2mA Ra 



2 ' 



0.3 + 0.91 



rriA \ 3 
GeVy 



2.28) 



When one can neglect both the form factor {Qa oo) and the relative velocity of the 
Sun with respect to the galactic halo {rja 0), Eq. ( p.27| ) simplifies to the form given in |0] 



Qa=oo 



KA{r,aji^^=o 



2 1 



1 — exp 



2 I ^esc(^) 



p. 



a 

~A 



(2.29) 
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III. GETTING SOME OF THE CAPTURED WIMPS OUT OF THE SUN 



The scattering events discussed in the previous Section create a population of solar- 
system bound WIMPs, moving (for a ~ 1 AU) on very elliptic orbits that traverse the 
Sun again and again . For the values of WIMP-nuclei cross sections we shall be mostly 
interested in below (corresponding to effective WIMP-proton cross sections (see section 
VII) in the range 4 x 10~^^ — 4 x 10~^^ cm^), the mean opacity of the Sun for orbits with 
small perihelion distances is in the range 10~^ — 10~^. This means that after only 10^ — 10^ 
orbits (i.e. ~ 10^ — 10^ yr) these WIMPs will undergo a second scattering event in the Sun, 
making them loose more energy, i.e. binding them even more to the Sun. From Eq. ( ^.11[ ) 
the average energy loss per scattering is 

(Q) = \^x(3t vL vUr) ^ (^) v^r) , (3.1) 

4 4 \mx / 

where we assumed rriA -C mx- 

Compared to the conserved energy before this (second) scattering event, E = mx a, 
the mean energy loss ( p.lD is quite large, because ('yesc('"))sun ~ (1000 km/s)^ ^ | a ~ | f |; ~ 
i (30 km/s)^. Even assuming a typical mass ratio as small as mig/l TeV ~ 1.6x 10^^ does not 
compensate the large factor 2{v1^^/v\ ~ 2 x 10^. The conclusion is that a second scattering 
event would typically bind the WIMP on an orbit of semi-major axis significantly smaller 
than 1 AU. This irreversible process {{Q) > 0) leads rather quickly (a few 10'^ — 10^ yr) to 
orbits of size comparable to the Sun (at which stage one might need to take into account the 
thermal velocities of the nuclei in the Sun, ultimately leading to a near thermal equilibrium 
between WIMPs and the core of the Sun |l^ ). 



The conclusion is that most of the population of WIMPs considered in the previous 
Section will end up quickly in the core of the Sun (where they will ultimately annihilate 
with each other, thereby generating an interesting indirect neutrino signal, see e.g. 0). The 
only way to save some of these WIMPs from this early demise is to consider the fraction 
of WIMPs that have perihelion distances rmin in a small range near the radius of the Sun 
Rs, say Rsi^ — e) < ''"min ^ Rs- As we argued in focussing on such a subpopulation 
of WIMPs has two advantages: (i) they traverse only a small fraction of the mass of the 
Sun and therefore their lifetime on such grazing orbits is greatly increased, and (ii) during 
this time, the gravitational perturbations due to the planets can build up and push them 
on orbits that no longer cross the Sun. We now tackle the latter perturbation problem. 

We first review in detail some concepts and notation of Hamiltonian dynamics. In stan- 
dard position-momenta variables the Hamiltonian describing the basic interaction between 
a WIMP and the Sun reads (r = |x|) 

H5(x,p) = ip2-f/(r), (3.2) 

where U{r) = +Gn J p(x') (i'^x'/|x — x'| is the (spherically symmetric) Newtonian potential 
generated by the mass distribution of the Sun. Thanks to the equivalence principle, the mass 
nix of the WIMP drops out of the problem (even when adding the perturbing influence of 
planets). Therefore in Eq. ( p.2|) and below we simplify the writing by factoring out mx 
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of all quantities, i.e. by working formally with a unit-mass WIMP. Because of the spher- 
ical symmetry of Eq. ( p.2|) we can first introduce spherical coordinates {r,9,ip,Pr,P9,p,p), 
with respect to which the problem is separable, and then work with some associated, con- 
venient action-angle variables ("Delaunay variables"). The Jacobi time- independent action 
SE{r, 9, ip) satisfying Hs{ 



i, dSo/dqi) 
SE{r,9,^) = Srir) + Sei9) 



E, is of the form 11 



= Jdr^2E + 2Uir)-^ + Jd9^J^-^^ + J,^. (3.3) 

Using Se we can introduce the usual (action-angle) Delaunay variables, traditionally denoted 

(L, G, H; i, g, h) The action variables L, G, H are related to E, J and appearing in 
Eq. ( p.3|) through (with pr = dSr/dr, pe = dS0/d9, p^ = dS^/dip) 

H = — p^d(p = J;,, (3.4a) 
Ztt Jo 

G = — pgd9 = J, (3.4b) 



2 /■'■max 2 /■'■max / Q2 

L = G+— Prdr = G + — drpE + 2U{r)- — , (3.4c) 



where the extra factors 2 compensate for the integrations on the intervals [6*1111115 ^max], 
['"min, '"max], Corresponding to half a period for these variables. The angle variables (with 
period 2tt) corresponding to L, G, H are respectively denoted h. [Their meaning will 
be discussed further below.] In these variables the Hamiltonian depends only on L and G 
(as is clear from Eq. (|3.4c| )), Tis = Ti-s{L, G), so that the general evolution equations. 



M_ m dg_ m cm_ m 

dt ~ ' 'dl ~ ^dG ' 'dt ~ ^dH ' ^ 

dL _ _dn dG _ _m dH_ _ _dH 

'dt ~ ~~di ' ~di ~ ~~d^ ' IT ~ ~~dh ' ^' ^ 



tell us, in the case of the problem (|3.2| ) that the action variables L, G and H are constant, 
while, among the angle variables, h is constant, but i and g evolve linearly in time: i = nt+ip, 
g = ut + go- Here, n = 2tt/P is the mean angular frequency of the radial motion (P = 
perihelion to perihelion period), and uj is the mean rate of advance of the perihelion. 

The central point is the following. A WIMP orbit with a generic perihelion distance 
'"mill ~ Rs will undergo a large perihelion precession Ac<j ~ 27r per orbit, i.e. cj ~ n, because 
the potential U{r) within the Sun is very modified compared to the exterior 1/r potential 
leading (by "accident") to the absence of perihelion motion. In other words, the trajectory 
of the WIMP will generically be a fast advancing rosette. This means that both angles ^ and 
g are fast variables. When adding in the small perturbing effect of the planets, i.e. when 
considering the total Hamiltonian, 

Htot = HsiL, G) + HpiL, G, H; i, g, h; Lp, . . . , ip, . . .) , (3.6) 
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where Tip (which contains a small factor jip = mpianct/^o) is the planetary perturbation, 
we can work out the (first-order) secular effects due to the planets by averaging over the 
fast variables £ and g (as well as the mean anomalies £p of the planets). Then the first 



two equations (|3.5b|) tell us immediately that the corresponding action variables L and G 
are secularly constant because planetary perturbations average out to zero (e.g. {dG/dt) = 
— {dH/dg)e^g = 0, because of the averaging over the fast angle g). As we shall see below L is 
essentially related to the semi-major axis a of the WIMP orbit, while G/L is related to the 
eccentricity e. The conclusion is that when the rosette motion is fast, planetary perturbations 
do not induce any secular evolution in the semi-major axis and in the eccentricity of the 
WIMP orbit. Therefore such WIMP orbits will necessarily traverse the Sun again and again 
and end up falling in its core. 

A different situation arises for WIMP orbits that graze the Sun. Throughout their orbits 
they feel essentially a 1/r potential due to the Sun, so that their rosette motion will be very 
slow. Consequently, the variable g will be slow for them (compared to i), and we cannot 



average on g. To tackle this problem we split the total Hamiltonian (|3.6| ) in three parts, 

T^tot = Ti-o + Hi + Hp , (3.7) 
where we take as unperturbed Hamiltonian the one corresponding to a point-like Sun, 



Ho = ^P^-^^, (3.8) 
2 r 



while the perturbations are 



Hi = -6U{r) ,np = -J2 Gn^p (- 

p \ 



Xx - Xp X, 



■Pi 



|3 



(3.9) 



Here, 6U{r) = U{r) — G^Mq/t is the non 1/r part of the potential generated by the Sun 
{6 U is zero when r > Rs and is responsible for the rosette effect when an orbit penetrates 
the Sun), and Hp denotes the planetary perturbations [|l3l. It is a sum over the planets 



with mass and heliocentric positions Xp. [The last term comes from the transformation 
between inertial (barycentric) coordinates and heliocentric ones.] We shall henceforth use the 
unperturbed Delaunay variables defined by the Hamiltonian Ho (corresponding to Keplerian 
motion). They are explicitly given in terms of the usual elliptic elements by 



'GNMQa , G= ^jGNMQa{l-e^) , H = ^G^v Mq a(l - e2) cosz, (3.10a) 
£ = mean anomaly, g = uj = periastron argument, 

h = VL = longitude of the ascending mode. (3.10b) 

Here, i denotes the inclination (with respect to the ecliptic), and we recall that the mean 
anomaly is, in Keplerian motion the angle £ = nt + io where n = 2tt/P is the radial circular 
frequency. It will be very convenient in this Section to use units such as 

1 



GnMq = 1^Ho = , L = , G= ^a(l -e2) , H = Gcosi. (3.11) 
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The fact that the unperturbed Hamiltonian depends only on the action variable L = Ir + Ie, 
Eq. (p.4c|) , is the famous degeneracy of the Coulomb problem. [In quantum language, L 
corresponds the principal quantum number Ug = + ng = rij. + j , while G = J corresponds 
to j and, H = Jz to m.] It implies immediately from the canonical equations (|3.5| ) that the 
only fast angle variable is i = nt + io+ perturbations, where n = dTio/dL = L^^ = a"^/^. 
This is just Kepler's law, = 1, in our units where Gn Mq = 1. 

We are interested in deriving the secular evolution of the elliptic elements a, e, i, or 
equivalently L, G, H, under the combined influence of the perturbations Hi and Hp. This 
is simply obtained by averaging the canonical equations over the fast angles, i.e. all the 
mean anomalies of the problem: i, ip. [We denote this average by an overbar.] By averaging 
Eqs. (|3.5b|) one sees easily that (in first order) L will be secularly constant (i.e. a = const). 



while G, H, g, h slowly evolve under the averaged perturbed Hamiltonian 

7^pert(^, G, H; g) = Hi{L, G) + Hp{L, G, H; g; Lp) . (3.12) 

As already indicated in Eq. ( p.l2| ), the averaged perturbed Hamiltonian depends (besides 
L = ^/a and Lp = ^Ja^ which can be treated as constants) only on G, H and g = uj. We 
consider, for simplicity, that the planets move on circular orbits. The lack of dependence 
on the last angle h = Q comes from the averaging over i and ip which establishes an 
azimuthal symmetry in the averaged Hamiltonian. Because of this azimuthal symmetry, 
dH/dt = —dHpert/dh = 0, i.e. H = Jz is a (secular) constant of the motion. Finally, we 
can treat L, H and Lp as constants, and consider only the evolution of the canonical pair 
{G,g) = {J,u!) under the Hamiltonian 7ipert(G, (?) = Hi{G) +Hp{G,g). But this is, in 
principle, a rather simple problem because we have now only one degree of freedom (one 
canonical pair), with a time- independent Hamiltonian. We conclude immediately that the 
perturbed energy is a constant of the motion. Therefore, it will be essentially enough to 
draw the phase-space picture of the problem via Hamiltonian level curves, 

HpertiG, g) = Hi{G) + HpiG, g) = E^,,, = const , (3.13) 

to be able to control the secular evolution of G, i.e. of the angular momentum J. We 
must now compute explicitly the values of Hi and Hp to see when planetary perturbations 
can, via the equation dG/dt = —dHpert/dg = —dHp/dg, force G = J to evolve ultimately 
towards large enough values of J, corresponding to orbits which no longer traverse the Sun. 

Let us first deal with the planetary perturbations. The calculation of Hp over the WIMP 
phase space is, in principle, a complicated matter because the WIMP can sometimes have a 
near collision with a planet, especially with Venus which is as massive as the Earth and nearer 
to the Sun. [Remember that we are interested in a population of WIMPs with a ~ 1 AU, 
so that they can be detected on Earth.] In fact, we tried to estimate separately the effect 
of such near collisions on the evolution of G = J. They will cause G to undergo a kind of 
random walk, which can certainly help to increase G beyond Js = Rs Vesc{Rs)- This effect is 
difficult to estimate with the kind of accuracy we attempt below for averaged perturbations, 
but it may be comparable to the effect we shall discuss below. By not including it explicitly, 
our estimate below of the fraction of WIMPs that can get kicked out of the Sun should thus 
be considered as a lower bound. To tackle analytically Hp we shall assume that the WIMP 
orbit is sufficiently smaller than the planetary orbits considered to be able to expand Hp in 
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powers of a/op and keep only the lowest significant (quadrupolar) contribution. [Actually, 
we have checked numerically that this quadrupolar approximation is surprisingly good even 
up to a ~ ap/2 for the very elliptic orbits we consider.] Doing the double average over £ and 
£p we find that the quadrupolar approximation to the second Eq. (p.9|) (actually the last, 
dipolar, term averages to zero) yields 



X,, 



E 



4a3 



(3.14) 



In terms of the eccentric anomaly u, we have |1T3 



u 



esinw , r = a (1 — ecosw) , z = sini (sina;5; + cosa;^) 
with {v denoting the true anomaly) 

X = r cos V = a (cos u — e) , y = r sin v = a -\/l — sin u . 



(3.15) 



(3.16) 



The average over i is easily computed as a modulated average over u, namely for any quantity 
q(i) = q(u) 



(9)^ = — (f dig = — (f du{l — e cosm) q = ((1 — e cosm) q)u ■ (3-17) 
2tt J 271 J 



This yields 



. 3 2 



sin^i [1 + (4 - Scos^cu)] , 



(3.18) 



so that 



i a 



3 2 3 



sin^ i 



■ sin^ i (5 sin^ cu — 1] 



(3.19) 



To use this Hamiltonian in the canonical equations ( |3.5| ) one should think of all the elliptic 
elements a, e, i, uj as being expressed in terms of L, G, H , g. In particular, we recall that 
L = y/a and H = G cos i can be considered as constants, while G = 



a 1 



and g 



Lu are 



the evolving variables. The dynamics resulting from the Hamiltonian Ti,p(G,g), considered 
separately from Hi{G,g), i.e. for orbits that stay outside the Sun, was studied some time 
ago by Kozai fl^ (see also the works of the Russian school in Ref. [^, and Refs. |1T6|JT7[] 
for studies in which this Hamiltonian is relevant). When considering the beginning of the 
evolution of G (while the WIMP is still in the outskirts of the Sun), it is essentially enough 
to use as Hamiltonian the (^-dependent part of Eq. (|3.19| ), i.e. 
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e sm % sm g . 



(3.20) 
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with e ~ 1 (very elliptic orbits) and i ~ const (because the fractional variation of the 
inclination while the WIMP is still in the Sun is small, being correlated through H = 
Gcosi = const to the fractional variation of G = J. Indeed, the WIMP just migrates 
through an outer skin of the Sun). 

Let us now consider the other perturbation, Tii, linked to the mass distribution in the 
Sun, Eq. (|3.9|) . For this term, it is more convenient to write the £-average in terms of an 
integral over the radial variable r (in terms of which 6U{r) is most easily expressed). To 
do that we can first replace the i (or time) average by an average over the area measure 
r"^ dv (Kepler's area law). Here, v is the true anomaly, so that the unperturbed trajectory 
on which we integrate is the ellipse r = a (1 — e^)/(l + ecosf). We can here replace a (1 — e^) 
by = (in our units) and replace e in the denominator (only) by one. [Indeed, we deal 
with very elliptical orbits with 1 — ~ -Rs/a <^ 1.] Differentiating the polar equation of 
the ellipse (approximated parabola) yields 



r dv 



Jrdr 



We then get 



^1 



vr a 



3/2 



r dr 



5f/(r) 



Rs f 



21/2 ^ ^3/2 J^^^^ .JT^ 



:6U{r) 



(3.21) 



(3.22) 



9 T ■ ~ 
I mm — 



Here, rinin(J) is the perihelion distance (closest radius) of the orbit such that 
('^ ^esc('"))^5 consistently with our definition of Tmin (for J = Jmm) in Section 2. [Indeed, 
■^escl^) — Gn Mq/t.] The maximum radius (which would be formally infinite for the 
parabolic approximation) plays no role because 6 U (r) vanishes by definition outside the 
Sun. 

To express explicitly Tii in terms of rmin (i-e. of G = J = -\/2r^) we need to know the 
radial dependence of 5 U{r) in the actual Sun. To do that we have fitted the radial mass 
distribution of the Sun, given numerically by Bahcall and Pinsonneault |jl8| , to simple power 
laws. We find that the mass distribution in the outskirts of the Sun (we are interested in 
the last outer few percent of the mass of the Sun) is well approximated by the fit 



6M{r) 



M, 







M(r) 




1 



(r > 0.55i?s) 



where 



0.0219 , = 0.907 i?c 



(3.23) 



(3.24) 



When using this fit, one must use Rs as an effective radius of the Sun beyond which the 
density vanishes (and 6U{r) = 0). [This neglects only the outer 0.13 % mass of the actual 
Sun.] From the (effective) radial mass distribution ( p.23| ) we can deduce both the density 
distribution (p(r) oc r~^Q{Rs — r)) and (by integration of dU{r)/dr = —GNM{r)/r'^) the 
Newtonian potential U{r) = Gn Mq/v + 6 U{r). This gives 



6U{r) 



GnMq 


'Rs 




(Rs^ 




Rs 


r 









e {Rs - r) 



(3.25) 
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By inserting (|3.25|) into ( p.22|) we finally get 



— e 'E^^^ r ■ .P 



where we used Gat Mq = 1, and where we defined the dimensionless function 



^(■^min) — O (1 ■^min) / 



lin 



1 



4 X 



(3.27) 



It is important to realize that as Xmin 1, i-e. as we get out of the Sun, the (Hamiltonian) 
function /i(xmin) tends to zero as fast as 

h{x^,n) ^1(1- X^nnf'^ © (1 " a^min) (3.28) 

5 

(because 5U{r) vanishes as (r — RsY)- By contrast, we shall see below that the efficiency 
of the outer skin in orbit-capturing WIMPs vanishes only as (1 — Xmin)^^^- This difference 
in asymptotic decrease will help in getting more WIMPs out of the Sun. 

As was explained above, the problem of selecting the values of the WIMP angular mo- 
mentum J = G for which planetary perturbations are strong enough to kick the WIMPs 
out of the Sun can be reduced to studying the level curves of the Hamiltonian (using only 
the crucial (/-dependent part of Hp, Eq. ( |3.20| )) 

2I, 



Ken = ^1 (G) + HliG, g) = ph[^]+a sin^ g , (3.29) 



where 



1 c; p^/^ 
a = -a smtJ2 — , (3=—^-^. (3.30) 



p 



Remember that the Hamiltonian function h{xj^in) ( p. 271) vanishes for Xmin > 1, i.e. when 
the angular momentum exceeds the value corresponding exactly to a grazing incidence, and 
behaves as (8/5) (1 — Xmin)^^^ when Xmin 1~- The level curves of ( p.29|) for Xmin < 1 



have the same shape as that of the Hamiltonian P'{G — GsY + asin^^f. This is nothing 
but a pendulum: impend = |Pe + \ ksiv?6. As is well known an essential element of the 
phase portrait of a pendulum is the separatrix, i.e. the level curve impend = | ^ which passes 
through the unstable equilibrium position pg = 0, 6 = 7r/2mod7r. This curve separates 
the oscillation motions (with 6 oscillating, nonharmonically, around ^ = mod tt) from the 
libration ones (with 6 changing monotonically, i.e. the pendulum going around in a sling 
motion). Similarly, in the case of the Hamiltonian ( |3.29| ) it is easy to see that the level curve 
Hpgrt = a (passing through G^ = 2Rs, g = 7r/2mod7r), i.e. 

/ r<2 \ 

Aa,i COSTS', (3.31) 




is a separatrix (for G < Gs = y2Rs). Here 
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Xa,^ = ^ = ^1 (^) si'^' ^ ' ai = 1 AU , (3.32a) 



The separatrix (|3.31|) divides between: (i) the uninteresting (for us) circulation motions 
where the "momentum" pe G — Gs keeps the same (negative) sign, i.e. the WIMP keeps 
traversing the Sun over and over again, and (ii) the oscillation (or libration) motions where 
the "momentum" pg G — Gs starts with a negative value and evolves so as to reach the 
value zero (corresponding to a trajectory which exactly grazes the Sun) in a finite time. 
Beyond the boundary G = Gs the Hamiltonian describing the evolution of the canonical 
pair {G,g) is the planetary perturbation Hp{L, H; G, g) given by Eq. ( p.l9| ) in which one 
must replace 

a = L'^ , = 1 - — , sin^ i = 1 - — , u = g , (3.33) 

where L and H are treated as constants. The corresponding phase portrait is drawn in 
Figure 1. [This Figure displays the level curves of the more exact perturbation Hamiltonian 
( p.l3| ), with Hp{L, H; G, g) given by Eq. ( |3.19| ), instead of the approximate expression ( p.20| ) 



used in the text.] The horizontal axis is the angle g (in radians), while the vertical axis 
is logio(G), where G = G/L = \/\ — denotes a "reduced" angular momentum. The 
numerical values used in Figure 1 are a = 0.844ai, and an initial inclination cos^(iin) = 1/3. 
The corresponding value for the Sun-grazing reduced angular momentum is Gs = 0.1. The 
separatrix discussed above lies in the domain G < Gs, i.e., logxo(G) < —1, and defines the 
dividing line (not explicitly shown) between the trajectories that stay always below Gs and 
those which evolve towards larger G values. Note that, starting on the "initial surface" 
G = Gsi most trajectories (except a small neighbourhood oi g = Omodvr) then evolve well 
away from the Sun to undergo large changes of G, up to values of order unity. In other 
words, once the very elliptic WIMP trajectory (initially with = 1 - ~ 2Rs/a ~ ^) 
exits the Sun, it undergoes, under the planetary perturbations, a slow, secular evolution 
of its eccentricity and inclination, up to values = 1 — ~ 1 and corresponding high 
inclinations « ~ f , keeping H/L = a/1 — cosz constant. 

We estimated the time scale for exiting the Sun, when starting on an initial "oscillation" 
trajectory with Gin < Gs- From the canonical equation dG/dt = —dHpcrt/^g and the 
constancy of H'-p^^^, Eq. p.29| ), the time it takes for Gin to evolve up to Gs is given by 
integrating 

^ = T («' - 4 [/3 hiG/Gs) - c]')-'/' , (3.34) 

where the constant c is related to the constant energy. By numerically integrating ( |3.34| ) 
over typical trajectories (not too near the separatrix on which it takes an infinite time to 
reach G = Gs) we found that it generically takes (for a ~ 1) less than 10^ WIMP radial 
periods (i.e. less than 10^ yr) for the eccentricity of the WIMP to increase sufficiently to exit 
the Sun. Then, when G > Gs the time scale for the evolution of G is given by the planetary 
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perturbations alone and is roughly [J2p /Up(a/ap)^]~^ longer than one orbital period, i.e. 
roughly 10^ yr for a ~ ai = 1 AU. After this time, the WIMP would, if it evolved only under 
the simplified planetary Hamiltonian Tip, come back again to low values of G, corresponding 
to Sun-penetrating orbits. Under the influence of Hi, it would then again bounce back 
away from the Sun in ~ 10^ orbits. For the scattering cross sections we shall be discussing 
below, the opacity of the small outer skin of the Sun we are interested in is typically smaller 
than ~ 10~^. Therefore the above process could persist for thousands of cycles before the 
WIMP gets scattered by a nucleus in the outskirts of the Sun. However, as we mentioned 
earlier, it is clear that the real gravitational interaction of the WIMP with planets is much 
more complicated than what is described by Tip. In particular, the non zero eccentricities 
of the planets, and the occurrence, once in a while, of a near collision with an inner planet 
will cause the elliptic elements of the WIMP to diffuse chaotically away from the simplified 
periodic history described above. Moreover, the very high eccentricities (for AU-size orbits) 
needed to traverse again the Sun represent only a very small fraction of the phase space into 
which the WIMP can diffuse. It thus seems clear that on time scales of several million years 
most of the population of WIMPs we are talking about will have irreversibly evolved onto 
trajectories on which the WIMP can survive (without being scattered again in the Sun) for 
the age of the solar system. We will come back later to the problem of the long-term survival 
of such WIMPs on orbits that stay within the inner solar system (rather than diffusing out 
into the outer solar system, and eventually to infinity). 



IV. ESTIMATING THE CAPTURE RATE OF LONG-SURVIVING, 
SOLAR-SYSTEM BOUND WIMPS 



In view of the previous estimates and arguments, we can consider that the population 
of WIMPs that, after a first scattering event in the Sun, diffuse out onto long-surviving 
solar-bound orbits is given by all the initial conditions Gin, 5'in, ^in which are "above" the 
separatrix (|3.31|) (meaning Ggeparatrix < G < G^). The logic for quantitatively estimating 
that population is the following. For each given initial values of a, and ii^, the separatrix 
defines, by solving Eq. (|3.31|) with respect to G, a corresponding minimum value of G 

Jmin — Gmin Gseparatrix('-^) fi'iii) ^in) ; ('^■■^) 

such that the trajectories with J > J^[^ end up out of the Sun. The differential capture rate 
corresponding to this class J > Jmin is then precisely defined by our previous result (p.22|) 
which defines 

= [Jmin(a,fi'in,iim),a] • (4.2) 

a a 

Then the actual capture rate is obtained by averaging ( [4.2|) over the distribution of initial 
values gin, ^in• 

Let us first approximate the result (p.22|) by a simpler expression. In the radial integration 
of Eq. ( p.22| ) the crucial features are the radial dependence of the abundance of element A, 
n/i(r), and the square root factor which vanishes at r = r^in- By contrast, the function 
KaIt, a), Eq. ( |2.27| ), varies fractionally very little over the small integration range Rs{l—e) < 
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r < Rs we are interested in. Therefore a good approximation of Eq. (|2.22| ) consists of taking 
out in front a factor KA{Rs,Oi) and performing the radial integration on the remaining r- 
dependent factors. To do that let us consider again the r-dependence of the total mass of 
the Sun. Denoting //(r) = M{r) / Mq, with < yu(r) < 1, we have from Eq. 



-R c dr dx 
dn{r)=3em — — = 3em — , (4.3) 



where was given in Eq. ( p.24| ) and where x = r/Rs- The distribution of the density of 
element A can be written as 

d x?T,yi(x) = 47rr drnA{r) = fA—^ dfi{r) , (4.4) 

where Ja denotes the mass fraction of element A in the outskirts of the Sun. Let us define 
the following ("capture") function 

/ \ dx r 3^ mill dx(^X XYnin) ^ / a r\ 

cix^in)= -Zi\^^--^= • (4-5) 

'-'-'iTkiii "-^ 111 ill 

In terms of this capture function, the capture rate reads 



dNA 



nx , Me 



, Vo ruA 

^min 



) , (4.6) 



dot 

where K\ is the "surface" value of the radial function explicitly defined in Eq. ( p.27| ) 

K\{a) = KA{Rs.OL). (4.7) 
Note that the surface value of the quantity Air) entering the error functions in ( p.27|) is 



/9- /_2 « 



{A{Rs)y = {\^aA)'^\%--^ \ ,n = GN Mq/Rs = (648.3 km/s)^ . (4.8) 



Actually, as said above, the a-dependence of is negligible (especially in view of all our 
other approximations) because a ~ f |; ~ (30km/s)^ so that v^/ct ~ 467 ^ 1. Therefore 
the a-dependence of dNA/da\j^_^.^_^ will come from the a-dependence of the "ejectable" radius 
a^min, to the estimate of which we now turn. 

Both the "capture function" c (xmin) and the previously introduced "hamiltonian func- 
tion" h (xmin), Eq. ( |3.27| ), can be explicitly expressed in terms of elementary functions. But 



these explicit expressions will not be really needed here. On the other hand, it is important 
to note the asymptotic behaviour of c {xmm) as Xjam ^ ■ 

C (a^min) ~ 2 •^min) ^ ■ (4'9) 

The fact that c(xmin) oc (1 — Xmm)^^^ decreases less fast than /i (xmm) oc (1 — Xmin)^^^, 
Eq. (p.28|) , as Xmin — * 1^ is important for us because the width of the separatrix ( p. 311) , 
^ (a^min) = 0{Xa,i), wiU be converted in a capture rate proportional to X^^^J, i.e. something 
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larger than the a priori expected small perturbation parameter Aa,j = a/ [3 (x Hp/a^- If we 
were to use only the asymptotic expressions ( F28D , (HD the width of the separatrix (^311) 



would be 1 — Xmin — (5 {a/aiY^^ (sini)''/^ (cos (7)^/^, and the corresponding capture 

rate would be proportional to 

Ca.yn.ptotic(x„.in) ^ \ (^-^^"' (^) ^ (sinz)^/^ {cosgf^ . (4.10) 
The actual capture rate will be larger than the one predicted by ( |4.10| ) because the actual 



function c (a;min) increases faster than (|4.9|) as one gets into the Sun. To combine some 
adequate numerical accuracy with the convenience of having analytical expressions we shall 
assume the approximate validity of the scalings in a/ai, sini and sin (7 predicted by Eq. (|4.10|) 



but calculate the precise numerical coefficient applicable for a/ai = sini = cosg = 1 by using 
the full numerical expressions of the functions h{x) and c(x), i.e. by inverting h{xi) = Ai 
in xi and computing c (xi). To do this we need the numerical value of Ai. First, taking into 
account the most significant planets, i.e. Venus, the Earth, Mars, Jupiter and Saturn, we 
find 

= 1.67xl0"^ (4.11) 

where we recall that ai denotes simply the basic unit for semi-major axes, namely the 
astronomical unit (AU). The other important numerical ingredients in Ai, Eq. ( p.32b ), are 




1/2 

^ ] = (236.9)^/2 ^ 15 39 ^ ^ o.0219 , (4.12) 




so that 



Ai = 0.0978. (4.13) 

The corresponding solution of h{xi) = Ai is xi = 0.729 (which means that we are 
typically dealing with the outer 27 % of the Sun radius-wise, containing in fact only ~ 2 % 
of the mass of the Sun). The corresponding value of the capture function is c(xi) = 0.172, 
so that the numerical combination effectively appearing in the capture rate ( [1.6D will equal 
(when a/ai = sini = cosg = 1) 

3 e„ c (a;i) = 0.0113. (4.14) 

This value must, according to Eq. ( |4.10| ), be scaled by (a/ai)^'^ (sini)^/^ (cos c/)^/^ Here, 
i and g are the initial values of the inclinations and perihelion argument. These quantities 
are random variables: g is expected to have a uniform distribution over [0,27r], while it is 
cosi which is expected to have a uniform distribution over [0, 1] (indeed, the direction of 
the vectorial angular momentum J is expected to be random on the celestial sphere). The 
averaging over these variables brings a factor 

((sini)^/5)cosi {{cos gf^^)g = 0.7567 x 0.6007 = 0.4545. (4.15) 
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Together with ( |4.14|) and the a-scahng we end up with a fraction kicked out of the Sun, 



(a) = 3 e„ (c ~ 0i (^^ j " , 0i ^ 5.13 x lO'^ . (4.16) 



Finally, if we define the A-dependent combination 



g^ = I±aAK%, (4.17) 



niA 



the rate (per a = Gn Mq/o) of solar capture of WIMPs that subsequently survive out of 
the Sun to stay within the inner solar system, reads 



dN. ' ■ ^ 



da 



a \ 10 nv 

- M^^gA. 4.18 



Note that the A-dependence is entirely contained in the quantity qa with dimensions 
[cross section] /[mass], e.g. [cm^]/[GeV/c^], or GeV~^ in particle units. The total capture 
rate is given by qa- 

A 

V. ESTIMATING THE PRESENT LOCAL PHASE-SPACE DISTRIBUTION OF 
SURVIVING SOLAR-CAPTURED WIMPS 



The last result ( [4.18|) of the previous Section gives the rate with which a fraction of the 
WIMP- Sun-scattering events populates a class of WIMPs that get out of the Sun with initial 
semi-major axis equal to a, and very high initial eccentricities such that 



l_e2 = 2^~8.44xl0-3^. (5.1) 

a a 

The question that remains is the following: what happens to this population while it 
slowly builds up during the 4.5 Gyr lifetime of the solar system? What is the present local 
distribution (in position and velocity space) as seen from the Earth of these WIMPs? These 
questions are very difficult to answer with precision because of the complexity of Hamiltonian 
dynamics in the solar system. One would need long-term numerical simulations to give 
reliable quantitative answers. However, we shall attempt here to make some estimates that 
will allow us to estimate the present observable effects of this population of WIMPs. 

There are two main worries about the long-term survival of this population. The first 
would be that they traverse the Sun again and again and end up getting accreted by it. 
We argued away this worry above ( because of the very small opacity of the outer skin of 
the Sun, of the repelling effect of the interaction with 6U{r), and of the small probability, 
given some additional chaos, that the WIMP again encounters the Sun). Note that existing 
asteroid simulations (see e.g. |jl9|) do not help in this respect because: (i) they restrict 
themselves to essentially planar initial data (while our WIMPs have fully three-dimensional 
trajectories), and (ii) they stop their numerical integrations as soon as an asteroid touches 
the Sun (while our WIMPs could survive ~ 10^ — 10^ passes in the outskirts of the Sun). 
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[Note that, in the case of the "lowest" (a = 2.1AU) orbits considered in |T^, 79% of them 



were ehminated because of impacting on the Sun.] A second worry concerns the possibihty 
that the adiabatic invariant a slowly evolves, in a quasi-random-walk, under the effect of 
near collisions with planets. This diffusion in a-space could lead a fraction of the WIMPs to 
have higher and higher values of a, possibly being ejected from the solar system. One can 
give a crude analytical estimate of the time scale on which a can change in the following 
way. Because of the exponential accuracy with which adiabatic invariants are conserved in 
absence of near collisions (i.e. in absence of near singularities in the complex plane, see e.g. 
||1 1|| ) , the cause of the random walk of a must be the existence of near collisions with some 
planet. Therefore, this effect will depend very much on the value of a. If a is smaller than 
aj/2 ~ 2.6 ai (where the label J stands for Jupiter) the WIMP orbit cannot cross Jupiter's 
orbit even if the eccentricity is very high. In that case, only one of the inner planets can 
have a near collision with a WIMP. Let generally be the mass of a planet whose orbit can 
cross that of a WIMP, with semi- major axis ax- The only small parameter in the problem of 
the non adiabatic evolution of ax is /ip = rrip/MQ (we use units such that G^Mq = 1). This 
non adiabatic evolution will be due to a more or less random succession of near collisions 
with the planet. Each collision would induce a velocity change 5f ~ fip/{bv) and an energy 
change Sax /ax ~ i/^pC-x/^ where b is the impact parameter. The rate of occurrence of 
such collisions is smaller as b decreases. Between two such quasi- collisions all the angular 
variables in the problem (which determine the ± sign in the energy change) have probably 
had the chance of being essentially randomized. After a long time t we can then consider 
that the total fractional change of ax, A a^, is a random walk so that one must consider 
(Aax/flx)^ ~ X! if^'pax/by. If we provisionally use units where ax = 1, one can see 

collisions 

that the typical time between two near collisions with impact parameter 6 is tf, ~ (6^ is 
an effective cross section around the planet, while the WIMP is evolving in a full 3D- volume 
~ a|^ = 1). The number of terms in the above random- walk is then iV ~ t/t;, ~ 6^ t so that 
(5ax)^ ~ N (fip/b)^ ~ fipt. Returning to ordinary units we find a typical diffusion law 

-±,tn = C^^;'Tx, (5.2) 
ax J to ^ 

where Tx is the orbital period corresponding to ax, and where C is some numerical constant 
of order unity. The dimensionless constant C is impossible to estimate with any accuracy 
on the basis of the previous rough argument (it can also contain a logarithm due to the 
integration over a relevant range of values of b). The estimate (|5.2|) suggests that the 
semi-major axis of Earth-crossing WIMPs diffuse on a very long time scale t/j ~ C (3.3 x 
10^)^ (ax/ai)^^^yr ~ C x 10^^ (ax/ai)'^/^ yr, so that we can essentially neglect the variation 
of ax over the age of the Sun. By contrast, the situation becomes dramatically different when 
c-x = aj/2 = 2.6 ai, because in this case it can cross the orbit of Jupiter with ^ip ~ (1047)"^. 
This leads to a much shorter diffusion time scale tj^ ~ 4.6 x C x 10^ (2 ax/oj)^''^ yr- The 
existence of such very different time scales depending on a^ < aj/2 or a^ > aj/2 is well 
known in asteroid research and is apparent in the results of long-term numerical simulations, 
see, e.g., Ref. [^. In principle such numerical simulations can give estimates of the diffusion 



times. It seems that a value of C ~ 0.1 is roughly compatible with several results and 
numbers in the literature PDI, pT]|, though the comparison might be difficult because asteroid 
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simulations start with quasi-planar initial conditions. 

To summarize: we expect, in first approximation, that the initial a-distribution derived 
in previous Sections can build up over ts = 4.5 Gyr with only small diffusion effects if 
a < aj/2 = 2.6 Oi, while if a > aj/2 this population is cut-off because of a fast diffusion in 
the outskirts of the solar system. As we said above, this conclusion, based on our analytical 
estimate Eq. ( |5.2|) , should be checked by dedicated long-term numerical simulations. 

Having discussed the secular evolution of a, we need now to discuss that of e and the 
other elliptic elements. As shown in Eq. (|5.1|) the initial values of the eccentricities are very 
near 1. The discussion of the phase portrait of the secular planetary Hamiltonian Tip in 
Section III showed that e and i undergo large oscillations with e evolving between values 
very near 1 and values of order one (say 0.3). The lower values of e depend on the value 
oi H = Jz = ^Ja^l — e^) cosi which is a secular invariant. Note that we now consider the 
dynamics implied by the Hamiltonian Hp, without the solar contribution Hi. The phase 
portrait of Hp differs from Fig. 1 in that the level curves escaping from the Sun and formerly 
librating around g = or tt, are now circulating, i.e. such that the angle g is monotonically 
evolving . The phase portrait of this pure Kozai Hamiltonian is represented, for instance, 
in ( Fig. 8), or (Fig. 3). For our case (a very small value of the constant of 
motion H/L) this means that most orbits will feature a small value of G/L, i.e. a large 
value of the eccentricity e, for a wide range of values of the periastron argument g ( a 
maximal eccentricity is reached for (7 = 7r/2 or 37r/2, but this maximum is broad, and large 
eccentricities are maintained during most of the evolution). 

In addition, the time-averaged probability for the eccentricity to fall in the range e± | de, 
will be peaked toward the extremal values of e(t) (because dt = de/e diverges there). The 
maximal value of e is very near 1 for all the WIMPs of the population, while the minimal 
value varies across the WIMP population, because it depends, among other initial data, 
on the value of the constant of motion H = a{l — e^) cosz. Therefore the overall time- 
averaged and population-averaged distribution function for e will have a peak only near 
e ~ 1. But, as we said above, this peak should be superposed onto a rather fiat distribution 
favoring high values of the eccentricity. ( Therefore the fact that the maximal eccentricities 
of each orbit are reached for g = 7i/2 or 37r/2, which implies a lack of spherical symmetry of 
the spatial distribution of the WIMPS, and which thereby somewhat disfavors the ecliptic 
plane, should not be numerically very significant). In the absence of detailed numerical 
simulations of the long-term evolution of our population of WIMPs this argument (based on 
the simplified secular Hamiltonian ( p.l9| )) suggests to use as an educated guess, for the mean 
distribution function of e, a distribution peaked at e = 1, i.e. simply a delta function S (1— e). 
Numerical simulations of asteroids, which go beyond the simplified Hamiltonian (|3.19|) and 
take into account near collisions with the planets, suggest that the analytically-expected 
large oscillations in e may actually be damped and lead to a population always having 
quite large eccentricities (see Fig. 4 in Ref. p^]). Resolving this question, and investigating 
the actual deviation from spherical symmetry of the WIMP population, would demand the 
running of long-term numerical simulations. In the present paper, we shall assume, for 
illustration purposes, a spherically symmetric population formally entirely concentrated at 
e ~ 1, which is technically simpler to deal with than a more realistic range of high values 
of e. Certainly, some of the details of the predictions that we shall sketch below depend on 
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this assumption, but we think that it is an appropriate approximation at this stage, though 
one which will have to be checked by dedicated numerical simulations. 

Under this assumption we can compute both the present space distribution and the 
present velocity distribution of our population of WIMPs. Let us first consider the spatial 
(numerical) density of WIMPs n{r). Consider first a subpopulation with some given values 
of a and e. By differentiating r = a (1 — e cos u), £ = u — e sin u we find 

di/dr = ±a-^r (e^ - (1 - r/af)-'/^ (5.3) 

so that the fraction of time, or elementary probability dp = 2\di/dr \ dr/2'K (where the extra 
factor 2 comes from the sign ambiguity), spent by this subpopulation within the radii r and 

r + dr is 

1 T dv 

dp= — - = (r - a (1 - e)) e (a (1 + e) - r) . (5.4) 

7ra2 ./e2_ (i_^/^)2 



From our previous arguments the number of WIMPs with semi-major axis within [a, a + 
da] is 

dN , dN , dN , 

—— da = —— da = ts —r- da , (5.5) 
da da da 



where ts — 4.5 Gyr is the age of the Sun, and where, from Eq. (^4.18 



^ = eQa,-a)0,(^) Me^,,o,, ,,o, = X:^?A. (5.6) 

The average number of WIMPs within the radii r, r + dr is obtained by multiplying ( p.4| ) 
and (|5.5| ) and integrating over a, 

dr N = A TT r"^ dr n{r) = J da dp, (5.7) 

so that the density of WIMPs reads (using \da\ = Gn Mq da/a^) 

1 , ^G„Me._diV9(.-a(l-e))e(.(l + e)-.) 

A-K^rJ a2 da a? ^e^-il-rjaf 

If we were to consider a population distributed in eccentricity with weight ^{e) de, we would 
need to add a further integration / de ip{e) in front of Eq. ( |5.8| ). Here, as we have mentioned, 
we shall simply assume that e ~ 1 for the entire population. It is then convenient to replace a 
by the new integration variable x = 2a/r. Inserting Eq. ( p.6|) in Eq. ( |5.^ ) and remembering 
the various step functions that limit the a-integration range we get 

1.9 



n{r) = u,nx[^) ' /„ (^) , (5.9) 



where 
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IrM^f (5.11) 

As /„ {aj/r) tends to a finite limit as r ^ aj, we see that the radial dependence of the 
WIMP density is essentially given by the factor (ai/r)^'^. This indicates that, if it were 
possible to do so, it would be easier to detect the WIMP population we are talking about 
nearer to the Sun, e.g. by building a detector in a mine (or examining WIMP-induced tracks 
in ancient mica) on Mercury! Let us also note that the radial distribution n{r) probably 
becomes cut off below some radius Vc < Qm = 0.387 ai, M being a label for Mercury. Indeed, 
for ax < O'm/'^ WIMPs do not have near collisions with any of the planets. Their secular 
orbital evolution should be rather well described by the quadrupolar Hamiltonian Tip, which 
means that they will episodically but repeatedly penetrate the outskirts of the Sun, thereby 
risking more to be scattered again. There should exist a critical semi-major axis a^, between 
Rs and aul'^-, such that when ax < cic the WIMP penetrates the Sun too often and finally 
gets accreted. 

In the following we shall focus on the value of the density of WIMPs at the orbital radius 
of the Earth, = ai = ^ AU. Let us define the enhancement in WIMP density due to the 
secondary population considered here as 

^ _ n{ai) _ (secondary) WIMP density at the Earth ^^-j 
nx halo WIMP density at infinity 



From Eqs. (|5.9|) , ( |5.10| ) one finds 

= 02 A 5(tot , (5.13) 
with (using Gn M^/ai = with ve = 29.78 km/s, and /„ (5.2) = 2.3474) 

02 = 01 = 01 X 0.0555 = 2.85 x 10"^ , (5.14) 

vl Mq 1.91x10^° ^ . 7.44x10^2 

A = -^ts^ = T-GeVcm-2 = — ^ geVf . 5.15 

Vo al i;o/220kms"^ t;o/220kms~^ ^ ^ ^ ^ 



Finally the local enhancement in density is 

-2 _ U-^-l-^ (-10) 



5.44 X 10^^ o 0.212 



= (../220kms-) ^ ^^"^ = (../220kms-) ' 

where g[ot^^ = 10^° (y'tot(GeV)^. The meaning of the 10^° (GeV)^ factor is that in ^ftot = 
{fA/rriA) ctaK^, with dimensions [mass]~^ x [cross section], one must express the mass 

A 

in units of GeV and the cross section in units of 10~^°GeV~^ (with h = c = 1). Note the 
conversion factor GeV-2 = 3.8938 x lO'^^ cm^ . We shall see in the next Section that g[ot'^^ 
can be higher than ~ 1, so that this new population could represent a significant increase 
above the standard halo WIMP density. 
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VI. OBSERVABLE SIGNALS FROM THE NEW WIMP POPULATION 



The secondary WIMP population discussed here could give rise to observationally signif- 
icant effects that have not traditionally been taken into account in the standard approach 
to dark matter detection, where one considers only the primary galactic halo population. 
The main observable signals from the new population are: (i) a new component, involving 
~ keV energy transfer, in the differential spectrum of direct detectors of WIMPs, (ii) a signif- 
icantly different angular spectrum in any detector with directional sensitivity, and and (iii) 
a possible significant increase in the indirect neutrino signal caused by WIMP annihilations 
in the Earth. To discuss the figures of merit associated with the new WIMP population, 
we need to fold in the velocity distribution of the WIMPs. Eq. (|5.3|) above, together with 
di/dt = n = {Gn Mq/o^Y^'^, shows that the radial velocity Vr = dr/dt of a WIMP passing 
at radius r reads 

v, = ±^-( ^E^) a' -{a- rff' . (6.1) 
r V a J 

The local velocity of a WIMP in the Earth frame is vioc = vx — v^;, where is the 
heliocentric WIMP velocity, whose radial component is ( |6.1|) . In our approximation where 
e ~ 1 for all the WIMPs \x is in the radial direction (with |vx| = \vt\)-, and therefore 
orthogonal to the Earth orbital velocity v^;. This yields ^i^^, = f ^ + f |; so that, with 
r = ai = 1 AU, 

v\oc = ve[^- • (6.2) 

With ai/2 < a < aj/2 = 2.6 ai, this predicts that the local velocity of the secondary 



WIMPs we discuss ranges only between ve = 29.8 km/s and — 1/2.Qve = 48.2km/s. 
These numbers depend on our approximation e ~ 1. They are, however, indicative of the 
values we might expect for the actual population. The local distribution is obtained by 
ehminating a using ( |5.5| ) and (|6.2|). 

Actually, the observable of most interest is the differential rate (events per kg per day 
and per keV) of scattering events in a laboratory sample made of element A 

FiiQ)r WQ)= L (6.3) 



Here, Q = q^/(2mA) = {m'^^^/mA) vf^^{l — cos^^cm) is the energy transfer from the WIMP 
X to the nucleus A, mj.ed{X,A) = uixmAl ijnx + rnA) is the reduced mass, n is the local 
number density of the considered WIMP population, F\ the form factor ( |2.7| ), v the local 
WIMP velocity and dn (f) the normalized speed distribution of the WIMP number density, 
with / dn{v) = 1. Note that for the standard WIMP population 7T,standard _ while for 
the new population discussed here n^^^ = Setix- Following Ref. 0, one can introduce for 
any WIMP population the dimensionless quantity 



TT 



T(e) = ^„„r (6.4) 
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For the standard galactic halo WIMPs 



T'standardlQ) ^ exp (-t^^in/^o) = ^Xp j_ "^-^^ j ^ (55) 

\ ^ "^o "^red / 

when neglecting the motion of the Earth with respect to the halo. For the low velocity 
WIMPs we are considering Tstandard(Q) — 1- Therefore the ratio 

- jdR/dQr- _ T-(Q) ^ 

= (^j:^/^g) standard = '^^ ystandard(g) " '^■^ ^ Wj • l^.b) 

This quantity is the figure of merit of most interest to us. It expresses the fractional increase, 
with respect to standard expectations, in the differential scattering rate. The distribution 



da 



da is obtained by taking the integrand of Eq. (|5.8| ) and normalizing it to one. Changing 



the integration variable from a to x = 2a/r = 2a/ai gives 



1 dxBix- 1)0(5.2 -x) 



One must then bring in the factor 1/v = 1/{ve V3 — 2x^^) from Eq. ( |6.2|) . Let us define 
the energy scale 

Qe^2^vI = 2 ( y m^vl. (6.8) 

For natural Germanium (m^) ~ 73 GeV this scale is Qe — 1-5 ijnx/ijnx + m^))^keV. Let 
us also define the function 

1 /■5-2 0(x — 1) 2 

= 77^ / . ^ 9 , T /X o -1 ' ^'"^'^^^^ = ^ 7, ' ^^-^^ 

i„(5.2j Jx^iM x^-^ \/x - 1 V3 - 2x ^ i- Q 

where /„(5.2) = 2.3474 is the integral ( |5.11| ). In terms of these definitions, the figure of 
merit ( |6.6|) reads 

p{Q)=P,d[^^^ , (6.10a) 

Pi = ^ ^5^ = 1.39 ^t^^ (6.10b) 
2 ve 

The function -D(g), with q = Q/QE is plotted in Figure 2. 

The plateau that is reached by D{q) as soon as g < 1 (i.e. Q < Qe) has value D{1) = 
0.803. This gives the maximal figure of merit [0 

piQE) = l.llg[;,''\ (6.11) 

If giot'^^ ~ 1 (see below), this yields a 100 % increase of the differential event rate below 
Q = Qij ~ keV. 
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Also, within our present rough approximation, e ~ 1, the direction of the incoming 
WIMPs from the new population will be strongly anisotropic. Indeed, not only are they 
entirely confined in the ecliptic plane, but even in this plane they have velocities whose local 
direction is within ±tan~"^ ~ 1/5.2) = ±51.8° of the vector — vg. This directional 
information might greatly help in distinguishing the real events from the background if 
one had a directional detector. Note, however, that long-term numerical simulations of the 
evolution of the elliptic elements of the WIMPs are probably needed to assess the robustness 
of this prediction, and that for the WIMP spectrum, based on the crude estimate e ~ 1. 
For instance, one expects the actual spectrum p{Q) to be a somewhat smoothed version of 
Figure 2, though the existence of a hump around Qe should survive. 

In view of this uncertainty on the exact spectrum of the WIMP population, we did not 
compute a precise figure of merit for the indirect neutrino signal from the Earth. Let us 
only point out the main features of the new signal. First, the capture by the Earth of a slow 
population v^^™ ~ ve, instead of the standard -y^^-ndard ^ jg j^iore effective. If we take only 
this effect in account, one would expect a figure of merit of order (vo/ve) ~ pi ~ 1.4 giot^^. 
However, another effect is also quite important. Because of the large ratio {vo/vf^^pj^ ~ 
(220 kms~^/ll kms~^)^ ~ 400 the Earth capture probability for incoming standard WIMPs 
is strongly peaked around the "resonances" mx — rn^ for some element A in the Earth . 
For the new population {vx /v^^^^^Y ~ 9 is much smaller, and the resonances become much 
broader. This means in particular that for masses mx > rn^Q (iron resonance) the neutrino 
signal will be much amplified, compared to standard expectations. However, the maximum 
WIMP mass for which capture is important depends very sensitively on the low-velocity 
cut-off Vc (measured in the Earth frame) of the WIMP population. Indeed, taking into 
account the fact that the escape velocity from the iron core of the Earth is v^^^ ~ 15kms~^, 
capture is only possible when > {vc/v^^^Y, which defines an upper bound on mx- For 
instance, within our present (rough) approximation, Vc = f © so that only WIMPs with mass 
mx < 2.62mFe — 147GeV would be captured by the Earth. 



VII. ESTIMATES FOR REALISTIC WIMPS 

To determine the relevance of the effects discussed here to the ongoing search for halo 
WIMPs, the actual numerical value of g[ot^^^ ^oi realistic WIMPs is of central importance to 
consider. To investigate this question we have explored the parameter space of the theoret- 
ically best motivated WIMP candidate: the lightest supersymmetric particle (assumed to 
be a neutralino) of the "Minimal Supersymmetric Standard Model" (MSSM) 0. Of course 
this is really not a single model, but a range of models, depending upon the assumptions 
one makes about such issues as Unification, and also the nature of Supersymmetry breaking. 
Because of this, detailed specific predictions of remnant neutralino densities, elastic scatter- 
ing cross sections, etc., are difficult to give with any generality. Indeed, our understanding 
of SUSY models is still developing, so that predictions of annihilation rates in the early 
universe, and thus remnant neutralino densities may require alteration [^]. 

In any case, for the purposes of this investigation it is worth exploring the general order of 
magnitude of predicted solar capture cross sections, and the resulting solar system density of 
SUSY WIMPs. To this end, to sample the many SUSY parameters, we have made use of the 
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specialized code Neutdriver written by Jungman, Kamionkowski and Griest [^]. This allows 
a calculation, using a specific parameterization of MSSMs, of annihilation rates, remnant 
neutralino densities, and elastic scattering cross sections on isotopes in the Sun, and in 
potential terrestrial detectors. 

Constraints on the SUSY parameter space are also model dependent, depending on 
whether one uses various GUT relations for gaugino masses, and also on assumptions about 
universality of scalar and fermion masses. These constraints are also evolving as new data 
from LEP, and from such processes as 6 — >■ s 7 are obtained. At the time the calculations 
reported here were performed, these constraints led us to sample the SUSY phase space 
described in Table 1 (conventions are those of Jungman etal. 0). 

As implied by the data in Table 1, a total of 9600 different sets of SUSY parameters were 
initially chosen. Among these some combinations, reported by Neutdriver, were unphysical 
or phenomenologically unacceptable for a variety of reasons. These were culled, and the 
remaining allowed configurations were utilized to determine remnant densities and scattering 
cross sections. 

This residual phase space has some characteristics that are important to distinguish 
here. First, we include neutralino masses as large as 400 GeV. These masses are larger than 
conventionally displayed in constraints by ongoing direct WIMP detection experiments. 
However as we are interested here primarily in knowing how broadly relevant our results 
might be, we wanted to explore as large a region of phase space as possible — independent 
of model builders' or experimentalists' preferences. This factor also played a role in our 
choice of WIMP cosmic densities to include in this analysis. From the broad phase space 
that survived the above cuts, we then selected those models that resulted in a remnant 
density in the range 0.025 < Vtxh? < 1. This range again is somewhat broader than is 
conventionally chosen for Q. However, given the upper limit ^^Baryon^^ < -026 from Big 
Bang Nucleosynthesis it remains possible that the non-baryonic abundance could be as 
small as the lower bound quoted above and still at least marginally exceed the baryon density 
in our own galaxy. This relaxed choice of density restriction, combined with the higher mass 
range we consider, implies that we allow models with somewhat higher cross sections on 
terrestrial targets than is usually displayed when SUSY constraint diagrams are displayed. 
(It is generally true, for example, that those models with the lowest remnant density today 
have the highest elastic scattering cross sections, for reasons made clear in the introduction 
to this paper.). In any case, the results quoted here are meant to be indicative of what 
one might expect for realistic WIMPs, and since SUSY model predictions are themselves 
evolving, the detailed model results quoted here should be taken as indicative of the general 
order of magnitude of one's expectations for SUSY WIMPs. Nevertheless, in order to explore 
how restricting the remnant neutralino density will affect the range of g^s expected for SUSY 
WIMPs, we also considered subset of the parameter space in which 0.1 < VLxh? < 1. 

Each neutralino has two possible modes of scattering on targets, in both the Sun and in 
terrestrial detectors. Because neutralinos are Majorana particles the non-relativistic limit of 
the scattering cross section generically involves a spin-dependent piece. In addition, exchange 
of scalar particles can produce a scalar, spin-independent piece. This latter term, if present, 
generally dominates for large nuclear targets, because in this case the WIMP can scatter 
coherently off of the entire nucleus with a coupling proportional to the atomic number of 
the nucleus A. Thus the cross section goes as A^. Moreover, the cross section generically 
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involves as factor the square of the WIMP-nucleus reduced mass, {mArrix/ {jnA + mx)Y, 
which, for heavy WIMPs, also increases as the square of the mass of the target nuclei. Since 
this latter quantity is also proportional to A, this implies that the scattering cross section 
for heavy nuclei can include a factor proportional to A^, which for nuclei as heavy as iron 
or germanium can be very large indeed. 

As a result, one generically finds that scattering cross sections are dominated by the 
coherent scalar piece in all cases except where this piece is suppressed due to various model- 
dependent factors. Moreover, with the exception of hydrogen and nitrogen, all other nuclei 
in the Sun are even-even nuclei, and for these nuclei the spin-dependent amplitude vanishes. 
In calculating the solar capture rate described in this paper, we included all known ele- 
ments in the Sun, with abundances given by current solar model calculations. All elemental 
abundances except for hydrogen and helium are taken from Jungman etal. H while the 



former two abundances are taken from Bahcall and Pinsonneault [0. We display the mass 
fractions used for the various elements in Table 2. 

In addition to calculating the capture factors qa, we also determined the scattering 
cross sections on germanium, which is currently the target material of choice in cryogenic 
detectors. In order to express the results in a more target independent way, however, we 
adopt a standard presentation of this cross section in terms of the effective WIMP-nucleon 
cross section. This is obtained by scaling from targets of atomic number A, using the 
assumption of coherent scattering, and is given by: 

aA, mx+mA , m^m, , ^ a a (m^ + tua? 
^ rnxmA mx + rap {nix + rnpY 

We present our results in Tables 3 and 4 and Figures 3-6. In the tables, in addition to 
the value of g[ot^^ , we also display the value of S'Hsliar ' ■6'Fe i Qq . In the figures we 
display WIMP-Nucleon effective cross sections as a function of mass for different models, 
where models with /i > (/i < 0) are displayed in odd (even) figures. The size of the model 
point gives the range of the value of g[ot^^ calculated for this model. The different figures 
refer to different cutoff values for the WIMP remnant cosmic density. In the figures we also 
show approximate limits obtained from direct detection experiments on o"p [23| under two 



assumed values for the local halo WIMP density (p = 0.3GeVcm~^ and p = O.lGeVcm"^). 

Several features should be clear from these results. First, gtot values in excess of unity 
are clearly possible, implying that realistic WIMPs to which the next generation of direct 
detectors will be sensitive should be expected to have a solar-system density in the region of 
the Earth that is significant. Second, we note that there is in general a monotonic relation 
between cr^ and gtot, with however, wide dispersion. Approximately one has 

gL''^ ~ ^p/(6 X lO-^^cm^) . (7.2) 

Also note that in this case, the dominant single contribution to gtot comes from scattering 
on iron in the Sun, although the net contribution to gtot from the combination of lighter 
elements is of the order of 40% of the total. For < 0, an additional possibility arises, at 
least for the lowest values of Qxh"^ and for low mass WIMPs. In this case, solar capture 
can be dominated by spin-dependent scattering off hydrogen so that the germanium cross 
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section, and hence the effective cxp can be several orders of magnitude smaller than in the 
case of WIMPs for which the dominant scattering on heavy nuclei is coherent. 

As one increases the lower limit on the remnant WIMP cosmic density today, several 
effects ensue. First, as expected, increasing this lower cutoff tends to decrease the mean 
value of Qtot- Also, for yU < the low mass low ap branch of WIMP phase space rapidly 
decreases in size, disappearing completely by the time the cutoff on Qh'^ exceeds 0.1. Thus, 
if the cosmic density exceeds this value, then a large solar system population is in one to 
one correspondence with WIMPs that should be directly detectable in the next generation 
of detectors. 

Finally, we have examined what would be the effect of changing the average velocity 
dispersion of halo WIMPs. We considered a range 180 < Vo < 270, which encompasses 
most estimates for this quantity for our galactic halo, and found that the results in all cases 
changed by less than 10% compared to those quoted above. It is interesting that while 
the change was too slight to be significant, the direction of the change was not monotonic, 
but depended upon the mass of the WIMP and the dominant target atom in the Sun. For 
heavy WIMPs whose dominant scattering was on heavy elements, increasing Vo increased 
the capture factor gtot- 



VIII. CONCLUSIONS 

The results presented here are quite encouraging, and motivate a consideration of detec- 
tion schemes that might probe down to keV energy deposits by WIMP scattering. If this 
is possible, the observation of a rise in the differential event rate for low energy events, of 
the form we describe here, could provide a very useful discriminant that could demonstrate 
that any claimed WIMP signal at higher energy is, in fact, due to halo WIMPs. While it 
is challenging to consider obtaining sensitivity to such low energy events (and more impor- 
tantly reducing the background of noise for such events), this may be less daunting than 
attempting to achieve directional sensitivity, which is the alternative discriminant that has 



been discussed P4| , p5| . Of course, if one had directional sensitivity, the signal we discuss here 
should be even easier to disentangle from backgrounds, as we expect it should be extremely 
anisotropic, as described earlier. 

Nevertheless, the analytical results presented here are in some sense still preliminary. 
While we expect the general quantitative features of this new WIMP population will be 
well approximated by the results presented here, full scale numerical simulations of the 
WIMP orbits under consideration will be necessary to confirm the details of our results. In 
particular, knowledge of the anisotropy of the distribution, as well as its energy spectrum will 
require such simulations, and the results presented here should be considered qualitative in 
these regards. In addition, such simulations, which incorporate the presence of the planets 
and allow close encounters, will be necessary to confirm that the WIMP population we 
focus on here is indeed long-lived in the solar system, and that its spatial distribution is 
not critically different from the simple spherically symmetric, high-eccentricity one we have 
assumed. 

One area that has not been investigated in detail here, and which certainly warrants 
further investigation, is the implications of this new distribution for indirect WIMP detection 
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via annihilations in the Earth. As we described in the text, it is quite hkely that this signal 
could be significantly enhanced, especially for heavy WIMPs, compared to that which is 
calculated for halo WIMP capture by the Earth. We expect, in fact, that new bounds on 
SUSY phase space may be possible on the basis of such considerations, compared to existing 
limits from underground neutrino detectors. Such an investigation is currently underway. 
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TABLES 



/X < U 


M2. oU-oUU 


mJ steps j 




Mi,M3 determmed by GUI relations 






ji: -oUU - -OU 


\L\j Steps j 




tanp. A -4U 


(4 steps) 




: /U -oUU UeV 


(o stepsj 




^squark^ 4 - 64 X 10^ GeV^ 


(4 steps) 


At > 


M2: 80-800 GeV 


(10 steps) 




Ml , M3 determined by GUT relations 






Ai: 150 - 800 GeV 


(10 steps) 




tan/3: 2 -40 


(4 steps) 




: 70 -500 GeV 


(3 steps) 




<uark: 4-64xl04GeV2 


(4 steps) 



TABLE I. SUSY Parameter Space Sampled in Estimates 



El(nncii( 


At;(niiic Alass Nuiiil)c-r 


Mmss Fraction 


H 


1 


0.7095 


He 


4 


0.2715 


C 


12 


3.87 xlO-3 


N 


14 


9.40 xlO-^ 





16 


8.55 xlO-3 


Ne 


20 


1.51 xlO-3 


Mg 


24 


7.39 xlO-"^ 


Si 


28 


8.13 xlO"'' 


S 


32 


4.65 xlO-^ 


Fe 


56 


1.46 xlO-3 



TABLE II. Elemental Mass Fractions Used in Solar Capture Estimates 
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mx(GeV) 


5tot 


m 


go 






Ai > 


















0.037 


384 


2.480 


0.008 


.410 


1.532 


1.38E-40 




0.038 


384 


2.239 


0.007 


.370 


1.381 


1.25E-40 




0.036 


384 


2.219 


0.007 


.367 


1.370 


1.24E-40 




0.032 


306 


1.943 


0.008 


.354 


1.136 


1.27E-40 




0.029 


306 


1.886 


0.007 


.342 


1.103 


1.24E-40 




0.027 


306 


1.871 


0.007 


.340 


1.094 


1.22E-40 




0.076 


353 


1.230 


0.004 


.210 


.746 


7.28E-41 




0.074 


353 


1.173 


0.004 


.201 


.710 


6.94E-41 


/i < 


















0.039 


357 


0.63 


0.002 


0.11 


0.38 


3.70E-41 




0.039 


357 


0.61 


0.002 


0.10 


0.37 


3.55E-41 




U.Uoo 


OO ( 


U.DU 


U.UUz 


U.iU 


U.oi 


O.0oil/-41 




0.025 


32 


0.59 


0.000 


0.00 


0.00 


8.17E-44 




0.037 


276 


0.59 


0.002 


0.11 


0.34 


4.16E-41 




0.035 


276 


0.58 


0.002 


0.11 


0.33 


4.08E-41 




0.034 


276 


0.58 


0.002 


0.11 


0.33 


4.08E-41 




0.029 


196 


0.55 


0.003 


0.12 


0.29 


5.12E-41 



TABLE III. Largest glJ ' values for 0.1 > ^xh^ > 0.025 
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mx(GeV) 


5tot 


m 


90 


9Fe 




Ai > 


















0.135 


397 


1.11 


0.003 


0.18 


0.69 


6.06E-41 




0.150 


80 


0.97 


0.013 


0.23 


0.47 


2.07E-40 




0.126 


80 


0.96 


0.013 


0.23 


0.47 


2.06E-40 




0.211 


80 


0.86 


0.012 


0.20 


0.42 


1.84E-40 




0.143 


397 


0.74 


0.002 


0.12 


0.46 


4.07E-41 




0.136 


397 


0.74 


0.002 


0.12 


0.46 


4.04E-41 




0.110 


316 


0.64 


0.002 


0.12 


0.38 


4.10E-41 




0.103 


316 


0.62 


0.002 


0.11 


0.37 


3.97E-41 




0.165 


359 


0.46 


0.002 


0.08 


0.28 


2.69E-41 




0.170 


359 


0.43 


0.001 


0.07 


0.26 


2.54E-41 




0.164 


359 


0.43 


0.001 


0.07 


0.26 


2.50E-41 




0.123 


396 


0.39 


0.001 


0.06 


0.24 


2.14E-41 




0.124 


278 


0.37 


0.002 


0.07 


0.21 


2.63E-41 


/i < 


















0.101 


361 


0.26 


0.001 


0.04 


0.16 


1.49E-41 




0.104 


361 


0.24 


0.001 


0.04 


0.15 


1.40E-41 




0.107 


199 


0.21 


0.001 


0.04 


0.11 


1.92E-41 




0.148 


402 


0.18 


0.001 


0.03 


0.11 


9.61E-42 




n 1 /I n 
U.i4y 


A no 


U.io 


U.UUi 


n HQ 
U.Uo 


nil 
U.ii 


n K QT? A o 

y.ooiij-4z 




0.130 


321 


0.18 


0.001 


0.03 


0.10 


l.lOE-41 




0.135 


321 


0.17 


0.001 


0.03 


0.10 


1.06E-41 




0.134 


321 


0.17 


0.001 


0.03 


0.10 


1.06E-41 




0.168 


362 


0.16 


0.001 


0.03 


0.10 


9.17E-42 




0.169 


240 


0.16 


0.001 


0.03 


0.09 


1.24E-41 



TABLE IV. Largest giJ ' values for 1.0 > O.xh'^ > 0.1 
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FIGURES 




FIG. 1. Level curves of the perturbation Hamiltoiiian describing the secular evolution of the 
canonical pair {G,g). Note the divide between trajectories that always stay within the Sun and 
those that get out. 
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q = Q/Qe 

FIG. 2. Shape of the differential rate ( per keV) of the additional scattering events caused by 
the WIMP population considered here. This figure represents the function D{q), where q = Q/Qe 
is the energy transfer in units of the characteristic energy scale Qe (in the keV range). 
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WIMP MASS (GeV) 

FIG. 3. gtot as a function of the effective 
WIMP nucleon cross section, ap, and WIMP 
Mass, for ^ > 0, and assuming Qxh^ > 0.025. 
Hatched curves represent experimental upper 
limits assuming p = CSGeVcm"^ (lower) and 
p = O.lGeVcm"^ (upper). 
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FIG. 4. Same as Figure 3 but with /x < 
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FIG. 5. Same as Figure 3, but assuming 
^xh^ > 0.1 



FIG. 6. Same as Figure 5, but with p < 
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